u = bfu(k).data;

dgdu = f_dgdu(u);
u_lim = 10^6;

epsilon1 = 10^(-2);

dfdu_sign = zeros(m*N,1);
dhdu = zeros(m*N,m*N);
label1 = zeros(m*N,m*N);
label2 = eye(m*N);
for i=1:m*N
    if abs(u(i)) >= epsilon1
        dfdu_sign(i) = sign(u(i));
        label1(i,i) = 1;
        dhdu(i,i) = sign(u(i));
        if abs(abs(u(i))-u_lim) <= epsilon1
            label2(i,i) = 0;
        end
    end
end

cvx_begin
    variables dfdu(m*N) lambda(n) residual(m*N) mu_kkt(m*N)
    % minimize sum(abs(residual))
    minimize square_pos(norm(residual, 2))
    subject to
        residual == dfdu + dgdu'*lambda + dhdu'*mu_kkt
        -1 <= dfdu <= 1 % deactivate constraint to produce same results as first method
        label1*dfdu == dfdu_sign
        label2*mu_kkt == zeros(m*N,1)
        mu_kkt >= 0 % not necessary for some reason?
cvx_end

residual_l2norm = norm(residual,2);